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Abstract 



We treat the calculation of gravitational radiation using the mixed 
timelike-null initial value formulation of general relativity. The determina- 
tion of an exterior radiative solution is based on boundary values on a timelike 
worldtube T and on characteristic data on an outgoing null cone emanating 
from an initial cross-section of T. We present the details of a 3-dimensional 
computational algorithm which evolves this initial data on a numerical grid, 
which is compactified to include future null infinity as finite grid points. A 
code implementing this algorithm is calibrated in the quasispherical regime. 
We consider the application of this procedure to the extraction of waveforms 
at infinity from an interior Cauchy evolution, which provides the boundary 
data on V. This is a first step towards Cauchy-characteristic matching in 
which the data flow at the boundary V is two-way, with the Cauchy and char- 
acteristic computations providing exact boundary values for each other. We 
describe strategies for implementing matching and show that for small target 
error it is much more computationally efficient than alternative methods. 
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I. INTRODUCTION 



We report here on an important step towards the ultimate goal of constructing numeri- 
cal relativity codes that calculate accurately in 3D the gravitational radiation at future null 
infinity. By "accurately" we mean (at least) second-order convergent to the true analytic 
solution of a well posed initial value problem. Thus our goal is to provide an accurate and 
unambiguous computational map from initial data to gravitational waveforms at infinity. Of 
course, uncertainties will always exist in the appropriate initial data for any realistic astro- 
physical system (e.g. in a binary neutron star system, the data for the metric components 
would not be uniquely determined by observations). But such a computational map enables 
focusing on the underlying physics in a rigorous way. 

Most relativity codes are second-order convergent, but because of boundary problems 
the convergence may not be to the true analytic solution of the intended physical problem. 
In order to explain this point, and to give the idea behind our method, we first briefly review 
some aspects of numerical relativity. The predominant work in numerical relativity is for 
the Cauchy "3 + 1" problem, in which spacetime is foliated into a sequence of spacelike 
hypersurfaces. These hypersurfaces are necessarily of finite size so, in the usual case where 
space is infinite, an outer boundary with an artificial boundary condition must be introduced. 
This is the first source of error because of artificial effects such as the reflection of outgoing 
waves by the boundary. Next, the gravitational radiation is estimated from its form inside 
the boundary by using perturbative methods, which ignore the nonlinear aspects of general 
relativity in the region outside the boundary. For these reasons the numerical estimate of 
gravitational radiation is not, in general, convergent to the true analytic value at future null 
infinity. The radiation properties of the Robinson- Trautman metric will be used to illustrate 
this effect. 

An alternative approach in numerical relativity uses the characteristic formalism, in 
which spacetime is foliated into a sequence of null cones emanating from a central geodesic. 
This approach has the advantage that the Einstein equations can be compactified [I[ so 
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that future null infinity is rigorously represented on a finite grid, and there is no artificial 
outer boundary condition. However, it suffers from the disadvantage that the coordinates are 
based on light rays, which can be focussed by a strong field to form caustics which complicate 
a numerical computation 0. Also, to date, the characteristic initial value problem has only 
been implemented numerically for special symmetries • 

Our ultimate goal is a 3D Cauchy-characteristic matching (CCM) code, which uses the 
positive features of the two methods while avoiding the problems. More precisely, the interior 
of a timelike worldtube T is evolved by a Cauchy method, and the exterior to future null 
infinity is evolved using a characteristic algorithm; boundary conditions at T are replaced by 
a two-way flow of information across T. In relativity, under the assumption of axisymmetry 
without rotation, there has been a feasibility study of CCM ||[10|; see also 0. CCM has been 



successfully implemented for nonlinear wave equations and demonstrated to be second-order 
convergent to the true analytic solution (which is not true in a pure Cauchy formulation 
with Sommerfeld outer boundary condition) fllf . 

While CCM has aesthetic advantages, it is important to ask whether it is an efficient 
approach. The question can be posed as follows. For a given target error s, what is the 
amount of computation required for CCM compared to that required for a pure Cauchy 
calculation? It will be shown that the ratio tends to as e — > 0, so that in the limit of high 
accuracy the effort is definitely worthwhile ||12|| . 

Our first step towards CCM is Cauchy-characteristic extraction (CCE) and we will 
present a partial implementation of CCE in this paper. The idea of CCE is to run a pure 
Cauchy evolution with an approximate outer boundary condition. A worldtube V is defined 
in the interior of the Cauchy domain, and the appropriate characteristic data is calculated 
on T; then characteristic algorithms are used to propagate this gravitational field to future 
null infinity ||13|| . CCE is simpler than CCM to implement numerically, because in CCE 
the data flow is one-way (Cauchy to characteristic) whereas in CCM the data flows in both 
directions. Note that the advantage of computational efficiency applies only to CCM and 
not to CCE. However, we will show that the advantage of second-order convergence to the 



true analytic solution does apply, under certain circumstances, to CCE. 

The work in this paper is part of the Binary Black Hole Grand Challenge, which is 
concerned with the gravitational radiation resulting from the in-spiral and coalescence of 
two arbitrary black holes. However, the methods described here are not limited to black 
hole coalescence and could be applied to gravitational radiation from any isolated system, 
either with or without matter. 

In Sec. U, we present a formalism for 3D characteristic numerical relativity in which 
the coordinates are based on null cones that emanate from a timelike worldtube T (Recall 
that existing codes are in 2D with null cones emanating from a timelike geodesic) |§. 
The characteristic Einstein equations are written as a sum of two parts: quasispherical 
(in a sense defined below) plus nonlinear. The discretization and compactification of the 
Einstein equations, with the nonlinear part ignored, is discussed in Sec. |T|. A computer 



code has been written and in Sec. |V| this code is tested on linearized solutions of the 
Einstein equations, and extraction is tested on the nonlinear Robinson- Trautman solutions. 
The Robinson- Trautman solutions are also used to investigate the error of perturbative 
methods in estimating the gravitational radiation at null infinity. Sec. [V] uses the formalism 
developed in Sec. [TJ to estimate the errors associated with the finite boundary in a pure 
Cauchy computation. This leads to the result concerning computational efficiency of CCM 
stated above. In the Conclusion we discuss the further steps needed for a full implementation 
of CCE, and also of CCM, and investigate under what circumstances CCE can provide 
second-order convergence to the true analytic solution at future null infinity. We finish with 
Appendices on the null cone version of gauge freedom and linear solutions of the Einstein 
equations, and on a stability analysis of our algorithm. 

II. CHARACTERISTIC EVOLUTION IN 3D 

This is the first step towards a 3D characteristic evolution algorithm for the fully nonlin- 
ear vacuum Einstein equations. Here we treat the quasispherical case, where effects which 
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are nonlinear in the asymmetry can be ignored. Thus the Schwarzschild metric is treated 
exactly in this formalism. However, rather than developing an algorithm for the linearized 
equations on a given Schwarzschild background, we will approach this problem in a mathe- 
matically different way. 

We adopt a metric based approach in which each component of Einstein's equation has 
(i) some quasispherical terms which survive in the case of spherical symmetry and (ii) other 
terms which are quadratic in the asymmetry, i.e terms of 0(X 2 ) where A measures deviation 
from spherical symmetry. We will treat the quasispherical terms to full nonlinear accuracy 
while discarding the quadratically asymmetric terms. For example, if were a scalar function 
we would make the approximation 

e^V + d e <j)d e( f) « e^V (1) 

Although this breakup is not unique, once made it serves two useful purposes. First, 
the resulting field equations are physically equivalent to the linearized Einstein equations 
in the quasispherical regime. (In the exterior vacuum region, the spherical background 
must of course be geometrically Schwarzschild but the quasispherical formalism maintains 
arbitrary gauge freedom in matching to an interior solution). Second, the resulting quasi- 
spherical evolution algorithm supplies a building block which can be readily expanded into 
a fully nonlinear algorithm by simply inserting the quadratically asymmetric terms in the 
full Einstein equations. 



A. The null cone formalism 

We use coordinates based upon a family of outgoing null hypersurfaces. We let u label 
these hypersurfaces, x A (A = 2, 3), be labels for the null rays and r be a surface area distance. 
In the resulting x a = (u, r, x A ) coordinates, the metric takes the Bondi-Sachs form p^JT5[1 



ds 2 = - f e 2/3 ^ - r 2 h AB U A U B ^j du 2 - 2e 2p dudr - 2r 2 h AB U B dudx A + r 2 h AB dx A dx B , (2) 
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where h AB hsc = &c an d det{h A B) = det(qAB) = Q, with g^s a unit sphere metric. Later, for 
purposes of including null infinity as a finite grid point, we introduce a compactified radial 
coordinate. 

A Schwarzschild geometry is given by the choice f3 = /3{u), V = e 2/3 (r -2m), U A = 
and Kab — Qab- To describe a linear perturbation, we would set h A B = Qab + ^Iab and 
would retain only terms in 'Jab which were of leading order in the linearization parameter 
A. Here we take a different approach. We express 



in terms of a complex dyad qA (satisfying q A qA = 0, q A qA = 2, q A = q AB qB, with q AB qBC = 
S A ). Then the dyad component J = hABQ A Q B /2 is related to the linearized metric by 
J — ^IabQ q B /2. In linearized theory, J would be a first order quantity. The 2-metric is 
uniquely determined by J, since the determinant condition implies that the remaining dyad 
component K = h ab q a q b /2 satisfies 1 = K 2 — J J. Refer to |T(| for further details, especially 
how to discretize the covariant derivatives and curvature scalar of a topologically spherical 
manifold using the 5 calculus. 

Because the 2-metric also specifies the null data for the characteristic initial value prob- 
lem, this role can be transferred to J. Terms in Einstein equations that depend upon J to 
higher than linear order are quadratically asymmetric. We do not explicitly introduce A as 
a linearization parameter but introduce it where convenient to indicate orders of approxi- 
mation. 



The Einstein equations G^ v = decompose into hypersurface equations, evolution equa- 
tions and conservation laws. In writing the field equations, we follow the formalism given 



Qab = - (qAqB + QaQb) 



(3) 



B. Quasispherical approximation 



in |710. We find: 



A, = ^rh AC h BD h A B, r hcD, 



(4) 
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(r 4 e- 2 ^ AB [/f), r = 2r 4 (r~ 2 p A ) r - r 2 h BC D c h AB ,r (5) 
2e- 2/3 V r = TZ- 2D A D A f3 - 2D A (3D A (3 + r - 2 e - 2l3 D A {r 4 U A ) >r - X -r\~^h AB U A U B , (6) 

where D A is the covariant derivative and TZ the curvature scalar of the 2- metric h A B- 

The quasispherical version of (Q) follows immediately from rewriting it as /3 r = A^, where 
ATg = rh Ac h BD h A B, r hcD,r/l-6 is quadratically asymmetric. This defines the quasispherical 
equation 

Thus in this approximation, /3 = H(u, x A ) + 0(\ 2 ). For a family of outgoing null cones which 
emanate from a nonsingular geodesic worldline, we could choose coordinate conditions so 
that H = 0. Similarly, in Minkowski space, we could set H = for null hypersurfaces which 
emanate from a non-accelerating spherical worldtube of constant radius. In a Schwarzschild 
spacetime, due to red shift effects, H need not vanish even on a spherically symmetric 
worldtube. Thus H represents some physical information as well as asymmetric gauge 
freedom in the choice of coordinates and choice of world tube. 

We wish to apply the same procedure to equations (|5]) and @. In doing so, it is useful 
to introduce the 0(A) tensor field 

Cab = \h° D ( V A h DB + V B h AD - V D h AB ) (8) 

which represents the difference between the connection D A and the unit sphere connection 
Va, e.g. (D A — V ' a)v b = —C AB vc- In solving for U A , we use the intermediate variable 

Qa = r 2 e- 2 ?h AB U B . (9) 

Then (||) reduces to the first order radial equations 

(r 2 Q A ) r = 2r 4 (r" 2 /^ - r 2 h BC D c h A B, r , (10) 
U A = r- 2 e 2 ?h AB Q B . (11) 
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We deal with these equations in terms of the spin- weighted fields U = U A q A and Q = Q A q A . 
To obtain quasispherical versions of these equations, we rewrite ( |I0| ) and (|Tl| 



as 



(r 2 Q), r = -r 2 q A q BC Vch AB , r + 2r\ A (r" 2 ^) r + N Q , (12) 
U = r~ 2 e w Q + N v , (13) 



where 



N Q = q A [r 2 h BC (CS A h DB , r + CE B h AD , r ) - r 2 (h BC - q BC ) V c h AB>r ] , (14) 
Nv = r~ 2 e 2 \ A {h AB - q AB ) Q B . (15) 

The quasispherical versions obtained by setting Nq = in ([TJ) and Ny = in (|U|) then 
take the form 

(r 2 Q) ir . = -r 2 (6j + m) jr + 2r 4 5 (r~ 2 p) p , (16) 
£/ r = r- 2 e 2 ^Q, (17) 

in terms of the spin-weighted differential operator 5. Since Q yT and U >r are asymmetric of 
0(A), we use the gauge freedom to ensure that Q and U are 0(A). 

Since V = r in Minkowski space, we set V = r + W in terms of a quasispherical variable 
W. Then (|6|) becomes 

W, r = \e W n - 1 - e^Qe 13 + ^r~ 2 (r 4 (9*7 + W)) r + N w , (18) 



where 



JV W = -ePV A [{h AB - q AB ) V B e?] - ^r 4 e^h AB U A U B . (19) 



We set N w = in (|18|) to obtain the quasispherical field equation for W iT . 
Next, by the same procedure, the evolution equations take the form 

2(rJ) ur - (r-V(rJ) )r ) = -r" 1 (r 2 W) r + 2r- 1 e p Q 2 e p - (r^W) r J + N J} (20) 

where 



Nj = ^f-[- ^C c AB V c e? - h AC C c DB (r*U D ) r - (h AC - q AC ) V B (r 2 U c ) r + \r^h Ac h BD U c r U D r 
- l -r 2 h AB ,.D c U c - r 2 U c D c h AB , r + r 2 h CD h AD>T (D C U B - D B U C )}. (21) 



The quasispherical evolution equation follows from (|20| ) by setting Nj = 0. 

The remaining independent equations are the conservation conditions. For a worldtube 
given by r = constant, these are given in terms of the Einstein tensor by 

CCy^r = 0, (22) 

where £ M is any vector field tangent to the worldtube. This expresses conservation of £- 



momentum flowing across the worldtube [ 13| . These equations simplify when the Bondi 



coordinates are adapted to the worldtube so that the angular coordinates x A are constant 
along the d u streamlines. Then U = on the worldtube and an independent set of conser- 
vation equations is given (in the quasispherical approximation) in terms of the Ricci tensor 
by 

R r u = r~ 2 W )U - 2r- 1 /3 u - \r^mW + \ (W + 317) p =0, (23) 
2q A R r A = 3 ((r- x W) ir - 4r- 1 /3 - 2/3 )U ) + 3 (J u - J r ) - r 2 U, ru = 0. (24) 

In the context of an extraction problem it is assumed that the interior solution satisfies 
the Einstein equations, and therefore that the conservation conditions are automatically 
satisfied on the extraction worldtube. 

The above equations define a quasispherical truncation of the vacuum Einstein equa- 
tions. Because these quasispherical equations retain some terms which are nonlinear in the 
asymmetry, their solutions are not necessarily linearized solutions in a Schwarzschild back- 
ground. However, in the perturbative limit off Schwarzschild, the linearized solutions to 
these truncated equations agree with the linearized solutions to the full Einstein equations. 
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III. DISCRETIZATION OF THE EQUATIONS 

In this section we describe a numerical implementation, based on second-order accurate 
finite differences, of the equations presented in Sec. |I[ We introduce a compactified radial 
coordinate x = r/(R + r), (with R being the extraction radius), labeling null rays by 
the real and imaginary parts of a stereographic coordinate £ = q + ip on the sphere, i.e. 
i A = (q,p). The radial coordinate is discretized as Xi = xq + {i — l)Ax for i = 1 . . . N x 
and Arc = (1 — Xq)/(N x — 1). Here Xq = 1/2 defines a world tube of constant surface area 
coordinate. The point x^ x = 1 lies at null infinity. The stereographic grid points are given 
by qj = jA and pk = kA for j, k = —N^ . . . N% and A = 1/N%. 

The fields J, /3, Q and W are represented by their values on this rectangular grid, e.g. 
Jijk = J( u n,Xi,qj,Pk)- However, for stability (see Appendix 0), the field U is represented 
by values at the points x i+ i = Xi + Ax/2 on a radially staggered grid (accordingly UQ k = 
U(u n ,x i+ i , qjiPk) )■ For the extraction problem, it is assumed that the values of the fields 
and the radial derivative of U are known at the boundary. In the following discussion, it is 
useful to note that asymptotic flatness implies that the fields f3(x), U(x), W(x) = r~ 2 W(x) 
and J(x) are smooth at x — 1, future null infinity X + . 



Hypersurface equation for Q 

In terms of the compactified radial variable x, the quasispherical field equation for Q 
reduces to 

2Q + x(l - x)Q >x = -x(l - x)(SJ + 3K) iX - 4d{3. (25) 

We write all derivatives in centered, second order accurate form and replace the value Qj-i 
by its average (Qi + Qi-2)/2. The resulting algorithm determines Qi in terms of values of J 
and (3 at the points Xi, and x;_2 

Qi + Qi-2 + Xi-x{\ - Xi - 1 ) ^ 1 2^~ 2 = -Xi-d 1 - x i-i) (® Jl 2Ax~ 2 + ^ 2Ax ~ 2 ^) ~ 4 ^ -1 ' 

(26) 
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(Here and in what follows, we make explicit only the discretization on the radial direction 
x, and we suppress the angular indices j, k). Since Eq. ( f26|) is a 3-point formula, it can not 
be applied at the second point, however, a suitable formula for X2 is given by 

Q t + Q R + x c (l - x c )®\^ = -x c (l - x c ) {® J -\^ + S^V^l " 29 (A + fa) , 

ox \ ox ox J 

(27) 

where the value of Qr is trivially obtained from the knowledge of U r at the boundary, and 
xc = {xi + xr)/2, 5x = Xi — xr. After a radial march, the local truncation error compounds 
to an 0(A 2 ) global error in Q at X + . 



Hypersurface equation for U 

In terms of the compactified radial variable x, the quasispherical field equation for U 
reduces to 

V f = -£. (28) 

We again rewrite all derivatives in centered, second order form. Because of the staggered 
placement of U, the resulting discretization is 

Ui = Ui-i + —j^^Ax. (29) 

The value of U at the first point is evaluated from the expansion 

Ui = U\ R + U, x \ R (x i+ , - x R ) + 0(A 2 ) (30) 

at the boundary. This leads to an algorithm for determining U at the point x i+ i in terms 
of values of Q at the points Xr lying on the same angular ray. After completing a radial 
march, local truncation error compounds to an 0(A 2 ) global error in U at X + . 
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Hypersurface equation for W 



The quasispherical field equation for W (|T%D, reexpressed in terms of x and W = W/r 2 , 

is 

Rx 2 W, x +2R^—W = -e w TZ - 1 - e^SBe^ + -Rx 2 (W + W) + R-^— (W + W) . 
1 — x 2 4 ' > x 1 — x v 

(31) 

Following the same procedure as in Eq. (|26|) we obtain 

-(1 - ( -e 2ft ^ + ie 2 ^- 1 ^^! - 2 - e ft SBe ft - e^SSe^-^ 

2 2 \2 2 / 

We obtain a startup version of the above with the substitutions x { _\ — > xc, Ax — > 5x, 
noting that at the boundary U, x is given. The above algorithm has a local error 0(A 3 )in 
each zone. In carrying out the radial march, this leads to 0(A 2 ) error at any given physical 
point in the uncompactified manifold. However, numerical analysis indicates an 0(A 2 log A) 
error at X + . 



Evolution equation for J 

In discretizing the evolution equation, we follow an approach that has proven successful 
in the axisymmetric case || and recast it in terms of the 2-dimensional wave operator 

□ (2 V = e- 2/3 [2^-(^, r ), r ] (33) 

corresponding to the line element 

da 2 = 2l (jl n v )dx IJ, dx v = e 2(i du[—du + 2dr], (34) 

where 1^ = is the normal to the outgoing null cones and n M is a null vector normal 
inwards to the spheres of constant r. Because the domain of dependence of da 2 contains the 
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domain of dependence induced in the (u, r) submanifold by the full space-time metric 
this approach does not lead to convergence problems. 

The quasispherical evolution equation (pOf) then reduces to 



e^U^(rJ) =H, (35) 

where 

H = -r-\r 2 W) tr + 2r- 1 e /3 S 2 e /3 - (r~ l W), r J. (36) 

Because all 2-dimensional wave operators are conformally flat, with conformal weight —2, we 
can apply to (^) a flat-space identity relating the values of rJ at the corners P, Q, R and S 
of a null parallelogram A with sides formed by incoming and outgoing radial characteristics. 
In terms of r J, this relation leads to an integral form of the evolution equation, 

(rJ) Q = (rJ) P + (rJ) s - (rJ) R + ± Jdu drH. (37) 

The corners of the null parallelogram cannot be chosen to lie exactly on the grid because 
the velocity of light in terms of the x coordinate is not constant. Numerical analysis and 
experimentation has shown [jnj that a stable algorithm results by placing this parallelogram 
so that the sides formed by incoming rays intersect adjacent w-hypersurfaces at equal but 
opposite x-displacement from the neighboring grid points. The elementary computational 
cell consists of the lattice points (n, i, k, I) and (n, % ± 1, k, I) on the "old" hypersurface and 
the points (n + 1, i, k,l), (n + 1, i — 1, k, I) and (n + 1, i — 2, k, I). 

The values of r J at the vertices of the parallelogram are approximated to second order 
accuracy by linear interpolations between nearest neighbor grid points on the same outgoing 
characteristic. Then, by approximating the integrand by its value at the center C of the 
parallelogram, we have 

(rJ) Q = (rJ) p + (rJ) s - {rJ) R + ^Au (r Q -r P + r s - r R ) He- (38) 



As a result, the discretized version of (^) is given by 
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(rJ)^ 1 = T ((r J)^ 1 , (r J)S , (rJ)? +1 , (rJ)? , (r J)?_J + —Aw (r Q - r P + r s - r R ) H c 

(39) 

where J 7 is a linear function of the (rJ)'s and angular indexes have been suppressed. Con- 
sequently, it is possible to move through the interior of the grid computing (r J)™ +1 by an 
explicit radial march using the fact that the value of rJ on the world tube is known. 

The above scheme is sufficient for second order accurate evolution in the interior of the 
radial domain. However, for startup purposes, special care must be taken to handle the 
second radial point. In determining (r J)™^ 1 the strategy (|38|) is easily modified so that just 
two radial points are needed on the n + 1 level; the parallelogram is placed so that P and 
Q lie precisely on (n + 1, 1, and (n + 1, 2, i, j) respectively. Note that the calculation of 
Tic poses no problems, since the values of W, U, and U >r are known on the worldtube and 
the value of W^ r on the worldtube can be calculated by ( |T8D . 

In order to apply this scheme globally we must also take into account technical problems 
concerning the order of accuracy for points near X + . For this purpose, it is convenient to 
renormalize (|39| ) by introducing the intermediate variable $ = (rJ)(l — x) = RxJ . This 
new variable has the desired feature of finite behavior at 1 + . With this substitution the 
evolution equation becomes 

$ Q = \x Q AuHc + (* P - \x P AuHc) + \^ ($5 + \x s MHo) - \^ (*r + \xrMU 

(40) 

where all the terms have finite asymptotic value. 

IV. TESTS 

Some of the fundamental issues underlying stability of the evolution algorithm are dis- 
cussed in Appendix [CJ We have carried out numerical experiments which confirm that the 
code is stable, subject to the CFL condition, in the perturbation regime where caustics and 
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horizons do not form. The first set of tests consist of evolving short wavelength initial null 
data, with all world tube data set to zero. In this case, the world tube effectively acts as 
a mirror to ingoing gravitational waves. The tests were run until all waves were reflected 
and radiated away to I + . In particular, data with | J |~ 10~ 6 was run from from u = to 
u = 40, corresponding to approximately 10 4 timesteps, at which time it was checked that 
the amplitude was decaying. 

In the second set of tests, we included short wavelength data with amplitude 10~ 4 for the 
boundary values of /3, J, U, Q and W on the world tube (with compact support in time) as 
well for the initial data for J (with compact support on the initial null hypersurface) . Again 
the code was run for approximately 4500 timesteps (from u = to u = 25), at which time 
all fields were decaying exponentially. This test reveals a favorably robust stability of the 
worldtube initial value problem, since in this case the world tube conservation conditions 
which guarantee that the exterior evolution be a vacuum Einstein solution were not imposed 
upon the worldtube data. 

We now present code tests for the accuracy of numerical solutions and their waveforms 
at infinity. The tests are based upon linearized solutions on a Minkowski background and 
linearized Robinson- Trautman solutions. These solutions provide testbeds for code calibra- 
tion as well as consistent worldtube boundary values for an external vacuum solution. In 
addition, we use numerical solutions of the nonlinear Robinson- Trautman equation to study 
the waveform errors introduced by the quasispherical approximation. 

A. Linearized solutions 

Appendices |A] and [FJ describe how to generate 3-dimensional linearized solutions on a 
Minkowski background in null cone coordinates and their gauge freedom. To calibrate the 
accuracy of the code, we choose a solution of ( |B^ ) and ( [B7D which represents an outgoing 
wave with angular momentum / = 6 of the form 

$ = (d z ) 6 ^, (41) 
15 



where d z is the ^-translation operator. The resulting solution is well behaved above the 
singular light cone u — 0. 

Convergence was checked in the linearized regime by choosing initial data of very small 
amplitude (| J |~ 10~ 9 ). We used the linearized solution fl4"l~|) to give data at u — 1, with 
the inner boundary at R — 1, and we compared the numerically evolved solution at u = 1.5. 
The computation was performed on grids of size N x equal 128, 192, 256 and 320, while 
keeping N x = 4A^. Convergence to second order was verified in the L±, L2 and norms. 



B. Robinson- Trautman solutions 



The Robinson- Trautman space-times |2(J contain a distorted black hole emitting purely 



outgoing radiation. The metric can be put in the Bondi form 

ds 2 = -{JC - ^)du 2 - 2Wdudr - 2rW A dudx A + r 2 q AB dx A dx B , (42) 

where K = W 2 [1 — L 2 (log W)] , L 2 is the angular momentum operator and W(u, x A ) satisfies 
the nonlinear equation 

12<9 u (logW) = W 2 L 2 /C. (43) 

The Schwarzschild solution (for a unit mass black hole) is obtained when W = 1. More 
generally, smooth initial data W(uo,x A ) evolves smoothly to form a Schwarzschild horizon. 
The linearized solutions to the Robinson- Trautman equation (|43| ) are obtained by setting 
W = 1 + <p an d dropping nonlinear terms in (f>: 

12d u (f) = L 2 (2 - L 2 )<f). (44) 

For a spherical harmonic perturbation <fi = A(u)Ye m this leads to the exponential decay 
A = A(0)e' ue{e+1){e2+e ' 2)/12 . 

These linearized solutions provide analytic worldtube data for our evolution code, along 
with the initial null data J = 0. We have used this as a check of code accuracy in the 
perturbative regime off Schwarzschild. With this data, the code should evolve J to be 
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globally zero to second order in grid size. Of particular importance for the extraction of 
waveforms, this should hold for the value of J at X + . We have carried out such a test with 
a small extraction radius (R = 3m) and a linearized solution of the form: 

W = 1 + A $l[(e~ 2u Y 22 + e~ Wu Y 33 )} (45) 

with A = 10~ 5 . The error norm 

| \£j\ | 2 = [ 1 du f dnj 2 (46) 



o 



was determined by integration over a sphere at X + with solid angle element dfl, and with an 
integration time of u\ = 2. The convergence rate to the true value was found to be 0(A 1,92 ). 
We have also obtained second-order accurate numerical solutions to the nonlinear 



Robinson- Trautman equation (44). See Ref. [16] for numerical details. This allows us 
to check the discrepancy between exact waveforms and waveforms obtained by regarding 
the whole spacetime in the quasispherical perturbative approximation. We have based this 
comparison on initial data in modes 

W\ u=0 = l + \&[Y lm ]. (47) 

In order to supply some physical perspective, the nonlinearity of the initial data is best 

1 /2 

measured in terms of e = ((M — Ms)/M) , where M is the initial mass of the system 
and Ms is the mass of the corresponding Schwarzschild background. (Here, Ms = 1). We 
also calculate the percentage of the initial mass which is radiated away during the entire 
course of our simulations. The Bondi news function determines the mass loss and it is 
also an appropriate physical quantity to invariantly describe radiative waveforms. In the 
coordinates adopted here, the news function is given by |21| 



N(u B ,x A ) = - W _1 3 2 W (48) 

where the Bondi time ub measured by observers at X + is related to u by dus/du = W. 

For various initial modes, we have calculated the news function for the numerical solution 
of the nonlinear R-T equation ( N™ ) and compared it with the news function of the linearized 
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solution (iVf). As expected, for small values of e they agree up to second order in e. Figure 
[I] graphs the time dependence of A e = iV™ — Nf (for a representative angle) for a system 
initially in a / = 2, m = 2 mode, which is the dominant gravitational radiation mode for a 
spiraling binary system. The figure illustrates that A e scales with e 2 . However, for larger 
e, corresponding to a total radiative mass loss greater than 2.5%, this is no longer the case 
and a noticeable discrepancy arises. For instance, as illustrated in Figure ^|, the difference 
between quadratically rescaling A e and its actual value is about 40% for a mass loss of 4%. 

Hence, this indicates that not only the first order perturbation treatment but also the sec- 
ond order treatment is grossly inaccurate in this regime. Serious discrepancies arise between 
N™ and Nf for ranges in which the mass loss is not extreme. In fact, A™ reveals an oscilla- 
tory behavior qualitatively quite different from the pure decaying mode of the perturbative 
solution, which has serious implications for the tidal acceleration which the radiation would 
produce in a distant gravitational wave antenna. As measured by the radiative component 
of the Weyl tensor ^4, the tidal acceleration is given by the Bondi-time derivative of the 
news function. In contrast to the monotonic decay of the perturbative solution, the actual 
behavior of ^4 exhibits damped oscillations. For a Y22 initial mode, Figure |3| shows the 
drastic difference between the numerically obtained ^4 and the corresponding ^ calculated 
with the perturbative solution. 

Similar nonlinear oscillations arise from other choices of initial data. Some partial ex- 
planation of this behavior might be possible using second order perturbation theory for the 
Robinson- Trautman equation |22| but the full behavior would require perturbation expan- 



sions far beyond practicality. 

V. COMPUTATIONAL EFFICIENCY OF CCM 

This section is concerned with the computational efficiency of a numerical calculation of 
gravitational radiation from an isolated system, such as binary black hole. By "computa- 
tional efficiency" we mean the amount of computation A (i.e. the number of floating point 
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operations) for a given target error e. We will show that the computational efficiency of the 
CCM algorithm is never significantly worse than that of a pure Cauchy algorithm; and that 
for high accuracy the CCM algorithm is always much more efficient. 

In CCM a "3 + 1" interior Cauchy evolution is matched to an exterior characteristic 
evolution at a worldtube r« = constant. A key feature is that the characteristic evolution 
can be rigorously compactified, so that the whole space-time to future null infinity may 
be represented on a finite grid. From a numerical point of view this means that the only 
error made in a calculation of the gravitational radiation at infinity is due to the finite 
discretization A; for second-order algorithms this is 0(A 2 ). The value of the matching 
radius r« is important and it will turn out that, for efficiency, it should be as small as 
possible. However, caustics may form if tm is too small. The smallest value of tm that 
avoids caustics is determined by the physics of the problem, and is not affected by either 
the discretization A or the numerical method. 

On the other hand, the standard approach is to make an estimate of the gravitational 
radiation solely from the data calculated in a pure Cauchy evolution. The simplest method 
would be to use the raw data, but that approach is too crude because it mixes gauge effects 
with the physics. Thus a substantial amount of work has gone into perturbative methods 
that factor out the gauge effects using multipole expansions and estimate the gravitational 



field at infinity from its behavior within the domain of the Cauchy computation |23T|?5f . We 
will call this method waveform extraction, or WE. While WE is a substantial improvement 
on the crude approach, it ignores the nonlinear terms in the Einstein equations. The resulting 
error will be estimated below. 

Both CCE and WE are extraction methods. That is, they use Cauchy data on a world- 
tube T to estimate gravitational waveforms at infinity, and they have no back effect on the 
Cauchy evolution. In both methods there is an error (which is difficult to estimate) due to 
the artificial Cauchy outer boundary condition. The difference between CCE and WE is in 
the treatment of the nonlinear terms between T and future null infinity and in the trunca- 
tion of the perturbative multipole expansion at some low order. WE ignores the nonlinear 
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terms, and this is an inherent limitation of a perturbative method. Even if it is possible to 
extend WE beyond linear order, there would necessarily be a cut-off at some finite order. 
The quasispherical implementation of CCE incorporates all multipole contributions but also 
ignores the nonlinear terms. However, it is in principle straightforward to incorporate these 
terms into the code. A full implementation of CCE would do so, and the nonlinear terms 
would be treated without error. 



A. Error estimate in WE 

We assume that a pure Cauchy evolution proceeds in a spatial domain of radius re, and 
the extraction is computed on a worldtube T of radius R, with R < tb- 
The evolution equation (^) may be written: 

{ r J),ur — quasispherical part + — Nj (49) 



with the nonlinear term Nj given by ( |21| ) (Actually, Nj also implicitly contains contributions 
from (J NQdr)/r 2 and J Nydr, and from the quasispherical approximation of terms in (p0|), 
but these effects are all of the same order as Nj). The order of magnitude of various terms 
can be expressed in terms of a function c(w, x A ) (whose time-derivative is the news function); 
note that c is not a small quantity. The expressions are: 

2 

J = O(^), h AB - q A B = O(^), h AB , r = O(^), = O(^), 

Q = O(-), U = 0(4), C% = O(-), W = 0(- 2 ). (50) 



These estimates are obtained by the radial integration of the field equations in Sec. [11 B 



assuming that the background geometry is Minkowskian and that the Bondi gauge conditions 
are satisfied. Should this not be the case then constants of order unity would be added to Q, 
U and W, and the effect of this would be to amend (^) by adding terms to the quasispherical 
part so that it represents wave propagation on a (fixed) non-flat background. However, the 
order of magnitude of terms in the nonlinear part would not be affected. Thus there is no 
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loss of generality, and a significant gain in simplicity and transparency, in performing the 
error analysis on a Minkowskian background. 

It is straightforward to confirm that the nonlinear correction to fl2"0|) involves terms of 
order 0(c 2 /r 3 ) or smaller. WE estimates the waveform at future null infinity from data at 
r = R. This could be made exact (modulo the error introduced by truncating the multipole 
expansion) if the nonlinear part of (|20|) were zero. Thus the error introduced by ignoring 
iVjis 

f°° c 2 c 2 
e(c«) = {c, u )exact - {c, u )we = / 0{—)dr = O(-). (51) 

JR r ft 

In the case of the collision of two black holes, with total mass M and with c = O(M), the 
error is 0(M 2 /R 2 ) and it is tempting to say that if extraction is performed at R = 10M 
then the expected error of the WE method is about 1%. This would be quite wrong because 
there is no reason for the constant factor in 0(M 2 /R 2 ) to be approximately 1. 



B. Computational efficiency 

A numerical calculation of the emission of gravitational radiation using a CCM algorithm 
is expected to be second-order convergent, so that after a fixed time interval the error 

e = 0(A 2 ) ~ k±A 2 , (52) 

where A is the discretization length and k\ is a constant. On the other hand, the same 
calculation using WE must allow for the error found in (|5T|) , and therefore after the same 
fixed time interval there will be an error of: 

e = 0(A 2 , R- 2 ) ~ max(k 2 A 2 , |§), (53) 

where k% and k% are constants. 

We now estimate the amount of computation required for a given desired accuracy. We 
make one important assumption: 
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• The computation involved in matching, and in waveform extraction, is an order of 
magnitude smaller than the computation involved in evolution, and is ignored. This 
is justified by the 2D nature of the extraction and matching processes as opposed to 
the 3D nature of evolution. 

For the sake of transparency we make some additional simplifying assumptions (otherwise 
some extra constants of order unity would appear in the formulas below but the qualitative 
conclusions would be unaffected): 

1. The amount of computation per grid point per time-step, a, is the same for the Cauchy 
and characteristic algorithms. 

2. The constants k±, ki in the equations above are approximately equal and will be written 
as k. 

3. In CCM, the numbers of Cauchy and characteristic grid-points are the same; thus the 
total number of grid points per time-step is 



3A 3 ' 



(54) 



4. In WE, the outer boundary tb is at y^2R; thus the total number of grid points per 
time-step is 

8lfl3 (55) 



3A 3 



It follows that the total amount of computation A required for the two methods is: 



8irR 3 a 



A C cm - 3 ^ , A WE - 4 ■ (56) 

Thus the method which requires the least amount of computation is determined by whether 
tm > R or tm < R- (Because of the assumptions (1) to (4) this criterion is not exact but 
only approximate.) 
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As stated earlier, the value of tm is determined by the physics, specifically by the condi- 
tion that the nonlinearities outside must be sufficiently weak so as not to induce caustics. 
The value of R is determined by the accuracy condition fl53|), and also by the condition that 
the nonlinearities outside R must be sufficiently weak for the existence of a perturbative 
expansion. Thus we never expect R to be significantly smaller than r^, and therefore the 
computational efficiency of a CCM algorithm is never expected to be significantly worse 
than that of a WE algorithm. 

If high accuracy is required, the need for computational efficiency always favors CCM. 
More precisely, for a given desired error e, Eq's. (|5^) and (|53D and assumption (2) imply 



Thus 



A = y^/k, R = ^hfe. (57) 



A-CCM — ^ > Si-WE — % £ 7/2 ' V / 



so that 



Accm = r 3 M e 3 / 2 
Awe kl 12 



as e 0. (59) 



This is the crucial result: the computational intensity of CCM relative to that of WE goes 
to zero as the desired error e goes to zero. 



VI. CONCLUSION 

The computer code described in this paper is a partial implementation of CCE. That is, 
given data on an r = constant worldtube T, the code calculates the gravitational radiation 
at future null infinity in the quasispherical approximation. A full implementation of CCE is 
currently being developed which addresses the following issues: 

• The ignored nonlinear terms in the Einstein equations must be calculated, discretized 
and incorporated into the code. 
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• Algorithms need to be developed to translate numerical Cauchy data near T into 
characteristic data on T. 

• In general F will be described in terms of Cauchy coordinates, and will not be exactly 
r = constant; the characteristic algorithm needs amendment to allow for this. 

Once a fully nonlinear CCE code has been achieved it will be possible, under certain 
circumstances, to obtain second-order convergence to the true analytic solution at future 
null infinity. For example, if T has radius R and the radius of the Cauchy domain is 
(> R), then causality implies that the gravitational field at T will not be contaminated by 
boundary errors until time tc ~ (re — R)> where t = at the start of the simulation. There 
is no analytic error in the characteristic computation, so there will be no analytic error in 
the gravitational radiation at future null infinity for the initial time period tc', under some 
circumstances this may be the time period that is physically interesting. Further, this time 
period tc may be extended by using results from the characteristic computation to provide 
the outer boundary condition in the Cauchy calculation. This would amount to a partial 
implementation of CCM since there would be data flow in both Cauchy to characteristic, and 
characteristic to Cauchy, directions (The implementation is only partial because R and r# 
are very different). Since the data flow is two-way, the possibility of a numerical instability 
arises. However, the timescale of the growth of any instability would be tc, and therefore 
such a computation could be safely run for a time of several tc] the results obtained would 
be second-order convergent to the true analytic solution. 

Once the technology for Cauchy to characteristic, and characteristic to Cauchy, data flow 
across an arbitrary worldtube has been developed, a full implementation of CCM will amount 
to taking the limit in which the outer boundary approaches the extraction worldtube. We 
are encouraged to believe that this is feasible, i.e. without numerical instability, because 



CCM has been achieved for the model problem of the nonlinear 3D scalar wave equation WL . 
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APPENDIX A: GAUGE FREEDOM 

Given a metric in a Bondi null coordinate system, the gauge freedom is 

5g ab = g ac d c e + g cb d c C - ?d c g ab (Al) 

subject to the conditions 5g 00 = 0, 5g 0A = and gAB$g AB = 0. These latter conditions 
imply the functional dependencies 

e = <u,x B ) (A2) 

£ A = f A (u, x B ) - J drg Ql g AB d B K (A3) 

and 

e = (r/2)(g 01 g 1A d A K-D B £ B ). (A4) 

For a spherically symmetric background metric we drop quadratically asymmetric terms to 
obtain 

f = 90 - e^r^dK (A5) 



and 



f 1 = - r - (§f + g£) = i [-r5§ (0 + 0) + 2e 2/3 g8ft] , (A6) 
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where qA^ A = £ and qAf A = 90, in terms of a complex scalar field (f>(u, x A ) 
This gives rise to the following gauge freedom in the metric quantities: 



and 



e 
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5 J = -3 2 + — S 2 k } (A7) 



r 



5e^ = - (e 2 ^) u - e 2 ^ r = - (A) u + Vs8(0 + 0), (A8) 



V e 2(3 

6u = e )U - -e,r - — $e m 



5V=-(2re + KV), u + Ve, r -re(^J . 



(A10) 



APPENDIX B: LINEAR SOLUTIONS 

We present a 3D generalization of a scheme [§] for generating linearized solutions off a 
Minkowski background in terms of spin-weight quantities a and Z, related to J and U 
by J = 25 2 a and U = dZ. We may in this approximation choose a gauge in which (3 = 
or otherwise use the gauge freedom to set (3 = H(u,x A ). In either case, W is given by the 
radial integration of the linearization of (|1^) and the remaining linearized equations reduce 
to 

(r 4 Z r ), r = -2r 2 (2 - L 2 )a r (Bl) 

and 

E:=2{ra) tUr -- r (r 2 a, r - l -r 2 Z^ =0, (B2) 

where L 2 = —35 is the angular momentum operator. 
Now set 

r 2 a r = (r 2 $) r (B3) 
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and 

r 2 Z jT = 2(L 2 - 2) $. (B4) 

Then 

E = rn$ + 2($ + a) )U -^(r 2 $) r + Z, (B5) 
where □ is the wave operator 

rU $ = 2 (r$) ur - (r$) „, + ^L 2 $. (B6) 



It follows that 



£,, = i(r 3 D$) r . (B7) 



Suppose now that $ is a complex solution of the wave equation □<& = 0. Then Eq. ( |B1| ) 
is satisfied as a result of ( |B~3| ) and (|B4|)and (|B~7| ) implies i? r = 0. If $ is smooth and 0{r 2 ) 
at the origin, this implies E = 0, so that the linearized equations are satisfied globally. The 
condition that $ = 0(r 2 ) eliminates fields with only monopole and dipole dependence so 
that it does not restrict the generality of the spin-weight 2 function J obtained. Any global, 
asymptotically flat linearized solution may be generated this way. 

Alternatively, given a wave solution $ with possible singularities inside some worldtube, 
say r = R, we may generate an exterior solution, corresponding to radiation produced by 
sources within the worldtube, by requiring E\r = or 



(B8) 



This is a constraint on the integration constants obtained in integrating ( |B3| ) and (Bl) which 
may be satisfied by taking Z\r — and 

a,„|*= f^(r 2 $) -O U- (B9) 



This determines an exterior solution in a gauge such that U\r = 0. 
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APPENDIX C: LINEAR STABILITY ANALYSIS 



In the characteristic formulation, the linearized equations form the principle part of 
the full system of Bondi equations. Therefore insight into the stability properties of the 
full evolution algorithm may be obtained at the linearized level. Here we sketch the von 
Neumann stability analysis of the algorithm for the linearized Bondi equations, generalizing 
a previous treatment given for the axisymmetric case. The analysis is based up freezing the 
explicit functions of r and stereographic coordinate £ that appear in the equations, so that 
it is only valid locally for grid sizes satisfying Ar << r and |A£| << 1. However, as is 
usually the case, the results are quite indicative of the stability of the actual global behavior 
of the code. 

Setting G = rJ and T = r 2 U and freezing the explicit factors of r and ( at r = R and 
( = 0, the linearization of the Bondi equations (0), ([H]) and (p0|) takes the form 



# 2 r rr -2T = -{RG tr -G) x (CI) 

and 

2 G <ur — G trr = — T^r,rC- (C2) 

a 

Writing C = s i + ^ s 2, introducing the Fourier modes G = e wu e lkr e tl isi^i2S2 ( w jtn real k, li 
and l 2 ) and setting Y = AG, these equations imply 

A — -i (1 - ikR) (h - il 2 ) I [2 (2 + k 2 R 2 )] (C3) 

and 

2w = ik - (l\ + ll) (1 - ikR) I [(AR) (2 + k 2 R 2 )} , (C4) 

representing damped quasinormal modes. 

Consider now the FDE obtained by putting G on the grid points rj and T on the staggered 
points r I+ i/2, while using the same stereographic grid (j and time grid u^. Let P, Q, R 
and S be the corner points of the null parallelogram algorithm, placed so that P and Q are 
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at level N + 1, R and S are at level N, and so that the line PR is centered about 77 and 
QS" is centered about r/ +1 . Then, using linear interpolation and centered derivatives and 
integrals, the null parallelogram algorithm for the frozen version of the linearized equations 



leads to the FDE's 
2 



R 

Ar ' I 1 '* 



r rj .3 - 2r /+ i 



+ r 7 . 



- r 



R 1 

^- {G1+1 - Gi) - - + Gi) 

(C5) 



(all at the same time level) and 

W+i ~ Lr i ~ W+i + °7 4Ar ^ _ /+1 ~~ °7-i ~~ ^7+2 + A ~ °7 ) 



Au 



r / pAT+1 pAT+1 , pAT 



1 r , 1 



(C6) 



where ^ represents a centered first derivative. Again setting T = AG and introducing the 
discretized Fourier modes G = e ™ ArA V fc/A V' lJlAsi e i ' 2j2As2 , we have 8 C = L and 5, 



where L = (l/2)[sin(Z 2 As 2 )/(As 2 ) + i sin(ZiAsi)/(Asi)], and 

2 



and 



reduce to 



A 



R 

Ar 



[1 — cos a) + cos a 



R /a\ 1 /a 

i-7— sin — cos — 

Ar \2J 2 V2 



(C7) 



and 



C — AD 



(C8) 



,C-AD 

where a = A;Ar, C = ze ia/2 sin(a/2) + (Au/4Ar)(l — cos a;) and D = {iLAu/AR) sin(a/2). 
The stability condition that Re(w) < then reduces to Re[C(AD — AD)] > 0. It is easy to 
check that this is automatically satisfied. 

As a result, local stability analysis places no constraints on the algorithm. It may 
seem surprising that no analogue of a Courant-Friedrichs-Levy (CFL) condition arises in 
this analysis. This can be understood in the following vein. The local structure of the 
code is implicit, since it involves 3 points at the upper time level. The stability of an 
implicit algorithm does not necessarily require a CFL condition. However, the algorithm is 
globally explicit in the way that evolution proceeds by an outward radial march from the 
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origin. It is this feature that necessitates a CFL condition in order to make the numerical 
and physical domains of dependence consistent. In practice the code is unstable unless 
the domain of dependence determined by the characteristics is contained in the numerical 
domain of dependence. It is important to note that if U (or Y) are not discretized on a 
staggered grid then the above analysis shows the resulting algorithm to be unconditionally 
unstable regardless of any CFL condition. 
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FIGURES 
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FIG. 1. A e for e\ = 0.14 and €2 = 0.22 (corresponding to a total mass loss of 0.6% and 1.2% 
respectively) for initial data in a Y22 mode. In this regime A e scales as e 2 , thus indicating that 
first order perturbation is valid in this regime. 
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FIG. 2. A e for ei = 0.14 and €2 = 0.4 (corresponding to a total mass loss of 0.6% and 4.6% 
respectively) for initial data in a Y22 mode. The difference between quadratically rescaling A e 
and its actual value is about 40%, indicating that second order perturbation is inaccurate in this 
regime. 

34 



1.0 




-2.0 1 1 1 1 1 1 1 1 1 1 1 

0.00 0.10 0.20 0.30 0.40 0.50 

u 

FIG. 3. \I>4 and ^ for a point lying 10 degrees above the equator and initial data in a Y22 
mode. The total mass loss is 4%. The insert shows the marked oscillatory behavior at early times. 
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